      SUBROUTINE PLANETS(TSEC50,KPLANET,KCOORD,KUNITS,POS,LUERR,IERR)
      IMPLICIT REAL*8(A-H,M,O-Z)                                
C
C PLANETS GIVES THE GEOCENTRIC POSITION OF A SPECIFIED PLANET AT TIME 
C TSEC50. POSITIONS ARE DETERMINED ANALYTICALLY USING KEPLERIAN MOTION
C ABOUT THE SUN, THEN TRANSLATING TO GEOCENTRIC.
C
C VARIABLE DIM I/O  TYPE  DESCRIPTION
C -------- --- ---  ----  -----------
C
C TSEC50    1   I    R*8  TIME AT WHICH THE POSITION IS WANTED.
C                         ACCURACY DOES NOT WARRANT ET-UT CONSIDERATION.
C                         IN SEC SINCE 1/1/50, 0.0 HR
C
C KPLANET   1   I    I*4  NUMBER IDENTIFYING THE PLANET.
C
C                          = 0, SUN
C                          = 1, MERCURY
C                          = 2, VENUS
C                          = 3, EARTH
C                          = 4, MARS
C                          = 5, JUPITER
C                          = 6, SATURN
C                          = 7, URANUS
C                          = 8, NEPTUNE
C                          = 9, PLUTO
C                          = OTHERWISE, ERROR. POSITION RETURNED AS ZERO
C
C KCOORD    1   I    I*4  COORDINATE SYSTEM TO WHICH POSITION VECTOR IS
C                         TO BE REFERENCED.
C
C                          = 1, GEOCENTRIC MEAN EARTH EQTR/EQNX OF 1950.
C                          = 2, GEOCENTRIC MEAN EARTH EQUATOR AND
C                               EQUINOX OF DATE(DATE=TSEC50)
C                          = 3, GEOCENTRIC TRUE EARTH EQUATOR AND
C                               EQUINOX OF DATE(DATE=TSEC50)
C                          = OTHERWISE, ERROR. POSITION RETURNED AS ZERO
C
C KUNITS    1   I    I*4  FLAG INDICATING THE LENGTH UNITS WANTED FOR
C                         POS.
C
C                          = 1, ASTRONOMICAL UNITS. 1 AU IS DETERMINED
C                               BY THE VALUE IN THE FUNCTION SUBROUTINE
C                               CONST.
C                          = OTHERWISE, KILOMETERS
C
C POS       3   O    R*8  CARTESIAN GEOCENTRIC POSITION VECTOR OF THE 
C                         PLANET CENTER.
C
C LUERR     1   I    I*4  FORTRAN UNIT NUMBER FOR ERROR MESSAGES.
C                          = ZERO OR NEGATIVE, NO MESSAGE IS POSSIBLE. 
C
C IERR      1   O    I*4  ERROR OUTPUT FLAG.
C                          =0, NO ERROR FOUND
C                          =1, ERROR IN INPUT. PLANET ID IS INVALID
C                          =2, ERROR IN INPUT. KCOORD VALUE IS INVALID
C
C**********************************************************************
C
C BY KAREN M. KEADLE, 8/83.
C     MODIFIED....    CJP, 12/84.....
C
C**********************************************************************
C 
      REAL*8 DEGRAD/57.29588051308232D0/
      REAL*8 TWOPI/6.283185307179586D0/
C
      PARAMETER NPLANETS = 9
      DIMENSION ELEMS(6,NPLANETS),PERIOD(NPLANETS)
C
      REAL*8 POS(3),HELIOPOS(3,2),KEPLER(6),DUM3(3)
C
      INTEGER INIT/0/
C
      IERR = 0
C
C FIRST-CALL LOGIC. LOAD ELEMS, THE KEPLERIAN ELEMENTS FOR THE PLANETS.
C COORDINATE SYSTEM IS HELIOCENTRIC, MEAN EARTH EQUATOR AND EQUINOX
C OF 1950.0. UNITS ARE KM AND RADIANS. EPOCH IS SEC SINCE 1/1/50,
C 0.0 HR. PERIOD IS IN SECONDS. SUNMU IS KM**3/SEC**2.
C
      IF(INIT.EQ.0) THEN
        INIT = 1
        CALL PLANETS1(EPOCH,ELEMS,PERIOD,SUNMU)
        END IF
C
C ERROR CHECKS
C
      IF(KPLANET.LT.0 .OR. KPLANET.GT.9) THEN
        IERR = 1
        POS(1) = 0.D0
        POS(2) = 0.D0
        POS(3) = 0.D0
        IF(LUERR.GT.0) WRITE(LUERR,9001) 'BAD PLANET NUMBER.',KPLANET
        GO TO 9999
        END IF
C
      IF(KCOORD.LT.1 .OR. KCOORD.GT.3) THEN
        IERR = 2
        POS(1) = 0.D0
        POS(2) = 0.D0
        POS(3) = 0.D0
        IF(LUERR.GT.0) WRITE(LUERR,9001)
     *      'BAD COORDINATE SYSTEM IDENTIFIER.',KCOORD
        GO TO 9999
        END IF
C
C
C SPECIAL CASE. EARTH IS ALWAYS AT 0,0,0.
C
      IF(KPLANET.EQ.3) THEN
        POS(1) = 0.D0
        POS(2) = 0.D0
        POS(3) = 0.D0
        GO TO 9999
        END IF
C
C FIND HELIOPOS(-,1), THE POSITION VECTOR OF THE PLANET WRT THE SUN, AND
C FIND HELIOPOS(-,2), THE POSITION VECTOR OF THE EARTH WRT THE SUN.
C BOTH VECTORS ARE IN MEAN OF 1950.0, EQUATORIAL COORDINATES. 
C
      DT = TSEC50 - EPOCH    ! ELAPSED TIME SINCE EPOCH
C
      DO 15 IPL=1,2
      IF(IPL.EQ.1) KPL = KPLANET
      IF(IPL.EQ.2) KPL = 3
      IF(KPL.EQ.0) THEN   ! IE, THIS IS THE SUN.
        HELIOPOS(1,1) = 0.D0
        HELIOPOS(2,1) = 0.D0
        HELIOPOS(3,1) = 0.D0
        GO TO 15
        END IF
C
C   MEAN ANOMALY CHANGE BETWEEN EPOCH AND TSEC50.
      DELMA = TWOPI * DT/PERIOD(KPL)
C
C   MEAN ANOMALY OF THE PLANET AT THE REQUIRED TIME
      MEANANOM = EQVANG( ELEMS(6,KPL) + DELMA )
C
C   SET KEPLER, THE KEPLERIAN ELEMEMTS AT THE REQ'D TIME
      KEPLER(1) = ELEMS(1,KPL)
      KEPLER(2) = ELEMS(2,KPL)
      KEPLER(3) = ELEMS(3,KPL)
      KEPLER(4) = ELEMS(4,KPL)
      KEPLER(5) = ELEMS(5,KPL)
      KEPLER(6) = MEANANOM
C
C   CONVERT KEPLER TO CARTESIAN
      CALL TOCART(SUNMU,KEPLER,0,HELIOPOS(1,IPL),DUM3,6,IER1)
      IF(IER1.NE.0) STOP 'PLANETS. CODING ERROR. STOP 1.'
C
   15 CONTINUE
C
C
C CALCULATE PLANET POSITION IN GEOCENTRIC EQUATORIAL, MEAN OF 1950.0
C COORDINATE SYSTEM.  NOTE THAT PLANET1 INITIALIZED THE ELEMENTS IN
C HELIOCENTRIC, MEAN EARTH EQUATOR AND EQUINOX OF 1950.0, SO ONLY
C TRANSLATION IS NEEDED.
C
      FACTOR = 1.D0
      IF(KUNITS.EQ.1) FACTOR = CONST(60)   ! AU VALUE
      DO I=1,3
        POS(I) = (HELIOPOS(I,1) - HELIOPOS(I,2)) / FACTOR
        END DO
C
C ROTATE FROM MEAN EQTR/EQNX OF 1950.0 TO OUTPUT SYSTEM, IF NECESSARY.
C
      IF(KCOORD.EQ.2) CALL VECM50MDT(TSEC50,1,POS,POS)
      IF(KCOORD.EQ.3) CALL VECM50TOD(TSEC50,1,POS,POS)
C
C
 9999 CONTINUE
      RETURN
 9001 FORMAT(/,' PLANETS ERROR. ',A,'  VALUE=',I7)
      END
